New Approach to Solving the Radiative Transport Equation

We have proposed a novel method for solving the linear radiative transport equation (RTE) in a three-dimensional macroscopically homogeneous medium. The method utilizes the concept of locally rotated reference frames and can be used with an arbitrary phase function of a random medium consisting of spherically-symmetric microscopic scatterers. The angular dependence of the specific intensity written in the spatial Fourier representation is obtained as an expansion into spherical functions defined in reference frames whose z-axes coincide with the direction of the Fourier vector k. Coefficients of this expansion are obtained by numerical diagonalization of several k-independent tridiagonal matrices whose elements depend only on the form of the phase function. The inverse Fourier transform is then computed analytically. This results in a closed-form expression for the RTE Green’s function in infinite space. Further, the plane-wave decomposition of the 3D Green’s function is obtained. It is shown that the modes in this decomposition are the evanescent plane waves. These modes can be used to construct the solution to the boundary value problem in the slab and half-space geometries. c © 2005 Optical Society of America OCIS codes: (170.3660) Light propagation in tissues; (030.5620) Radiative transfer Image reconstruction in optical diffusion tomography (ODT) requires formulating a mathematical model that governs propagation of near-IR light in biological tissues [1–3]. A description based on Maxwell equations is too detailed and is not currently used. A more practical option is the use of the radiative transport equation (RTE) or the diffusion equation (DE), which is an approximation to the former. Currently, the vast majority of ODT implementations rely on the DE which is much simpler mathematically. However, the diffusion approximation can not be used in many practically important cases. In particular, it is not accurate in regions with relatively low scattering or high absorption (such as, for example, voids filled with clear fluids), near sources and boundaries, and in optically thin samples. If one of the above situations is encountered, propagation of the near-IR light in tissues must be described by the RTE. Unfortunately, the RTE is notoriously difficult to solve, even in homogeneous medium with simple boundaries. This is especially true in the case of highly forward-peaked scattering which is typically encountered in biological tissues. Numerical approaches of current importance for obtaining forward solutions to the RTE are based on the discrete ordinate method [4], including the Fokker-Planck approximation for sharply forwardpeaked scattering [5], cumulant expansion [6,7], modifications of Ambarzumian’s method [8,9], and different levels of the PL approximation [10, 11]. Solving the inverse problem for the RTE is an even more daunting task. Historically, the first such solution was obtained with use of Monte-Carlo simulations [12]. However, the Monte-Carlo process requires an extremely large number of random walkers to obtain statistically reliable solutions. A more commonly used alternative is the discrete-ordinate method (approximation of the angular integral in the RTE by a Riemann sum). This method was used, for example, in [13] in conjunction with an algebraic iterative technique for inverting the two-dimensional RTE. The three-dimensional inverse problem was considered in [14, 15]. A hybrid method in which the medium is assumed to be separated into highly scattering regions where the diffusion approximation is valid and non-scattering regions where ballistic (geometrical) propagation is valid was developed in [16]. This approach does not require direct solution of the RTE. A very promising emerging approach is based on the cumulant expansion [17]. An alternative to the discrete-ordinate method is the method of spherical harmonics, often referred to as the PL approximation in cases with special symmetry. This approach has the advantage of expressing the angular dependence of the specific intensity in a basis of analytical functions rather than in the completely local basis of discrete ordinates. However, when no special symmetry is present in the problem, the method a102_1.pdf MH2.pdf © 2006 OSA/BOSD, AOIMP, TLA 2006